TetMesh Geodesics#

[1]:
import numpy as np
from lapy import TetMesh
from lapy.plot import plot_tet_mesh
import plotly
# plotly.offline.init_notebook_mode(connected=True)
import plotly.io as pio
pio.renderers.default = "sphinx_gallery"
[2]:
T = TetMesh.read_vtk('../data/cubeTetra.vtk')
#T.is_oriented()
T.orient_()

--> VTK format         ...
 --> DONE ( V: 9261 , T: 48000 )

Flipped 24000 tetrahedra
[2]:
24000
[3]:
from lapy import Solver
fem = Solver(T,lump=True)

evals, evec = fem.eigs(10)

TetMesh with regular Laplace
Solver: spsolve (LU decomposition) ...
[4]:
# also get A,B (lumped), and inverse of B (easy as it is diagonal)
A, B = fem.stiffness, fem.mass
Bi = B.copy()
Bi.data **= -1
[5]:
evnum=1
cutting = ['x<0.5']
# also here we comment all plots to reduce file size
# uncomment and take a look
plot_tet_mesh(T,vfunc=evals[evnum]*evec[:,evnum],plot_edges=None,plot_levels=False,cutting=cutting,edge_color='rgb(50,50,50)',html_output=False,flatshading=True)

Mesh is oriented, nothing to do
Found 3040 triangles on boundary.
[6]:
from lapy.diffgeo import compute_gradient
from lapy.diffgeo import compute_divergence
grad = compute_gradient(T,evec[:,evnum])
divx = -compute_divergence(T,grad)
vfunc = Bi*divx

[7]:
cutting = ['x<0.5']
plot_tet_mesh(T,vfunc=vfunc,plot_edges=None,plot_levels=False,cutting=cutting,edge_color='rgb(50,50,50)',html_output=False,flatshading=True)

Mesh is oriented, nothing to do
Found 3040 triangles on boundary.
[8]:
np.max(np.abs(vfunc-(evals[evnum]*evec[:,evnum])))
dd = np.abs(vfunc-(evals[evnum]*evec[:,evnum]))
max(dd)
[8]:
0.006755829
[9]:
from lapy import heat

tria = T.boundary_tria()
bvert = np.unique(tria.t)

u = heat.diffusion(T,bvert,m=1)
cutting = ['x<0.5']
plot_tet_mesh(T,vfunc=u,plot_edges=None,plot_levels=True,cutting=cutting,edge_color='rgb(50,50,50)',html_output=False,flatshading=True)
Found 4800 triangles on boundary.
TetMesh with regular Laplace
Matrix Format now:  csc
Solver: spsolve (LU decomposition) ...
Mesh is oriented, nothing to do
Found 3040 triangles on boundary.
[10]:
# compute gradient of heat diffusion, normalize it, and compute the divergence of normalized gradient
tfunc = compute_gradient(T, u)

# flip and normalize it
X = -tfunc / np.sqrt((tfunc ** 2).sum(1))[:, np.newaxis]
X = np.nan_to_num(X)

# compute divergence
divx = compute_divergence(T, X)
[11]:
# compute distance
from scipy.sparse.linalg import splu
useCholmod = True
try:
    from sksparse.cholmod import cholesky
except ImportError:
    useCholmod = False

A, B = fem.stiffness, fem.mass # computed above when creating Solver

H=A
b0=-divx

# solve H x = b0
print("Matrix Format now: "+H.getformat())
if useCholmod:
    print("Solver: cholesky decomp - performance optimal ...")
    chol = cholesky(H)
    x = chol(b0)
else:
    print("Solver: spsolve (LU decomp) - performance not optimal ...")
    lu = splu(H)
    x = lu.solve(b0)

x = x - np.min(x)
Matrix Format now: csc
Solver: spsolve (LU decomp) - performance not optimal ...
[12]:
cutting = ['x<0.5']
plot_tet_mesh(T,vfunc=x,plot_edges=None,plot_levels=True,cutting=cutting,edge_color='rgb(50,50,50)',html_output=False,flatshading=True)
max(x), 0.5*np.sqrt(3.0)
Mesh is oriented, nothing to do
Found 3040 triangles on boundary.
[12]:
(0.69931483, 0.8660254037844386)
[13]:
#debug gradient
v1func =  T.v[:,0]* T.v[:,0] + T.v[:,1]* T.v[:,1] + T.v[:,2]* T.v[:,2]# x coord

grad = compute_gradient(T,v1func)
glength = np.sqrt(np.sum(grad * grad, axis=1))
fcols=glength
fcols=grad[:,2]
cutting = ['x<0.5']
cutting = None
plot_tet_mesh(T,vfunc=None,tfunc=fcols,plot_edges=None,plot_levels=False,cutting=cutting,edge_color='rgb(50,50,50)',html_output=False)
Mesh is oriented, nothing to do
Found 4800 triangles on boundary.
[14]:
divx = compute_divergence(T, grad)
divx2 = Bi * divx
cutting = ['z<0.5']
plot_tet_mesh(T,vfunc=divx2,plot_edges=True,plot_levels=False,cutting=cutting,edge_color='rgb(50,50,50)',html_output=False,flatshading=True)
Mesh is oriented, nothing to do
Found 3040 triangles on boundary.
[15]:
divx2[5000:5010]
[15]:
array([5.9999948, 6.0000215, 6.0000215, 5.999988 , 6.000053 , 5.999975 ,
       5.9999676, 6.000024 , 6.000013 , 6.000008 ], dtype=float32)
[16]:
T.avg_edge_length()
[16]:
0.06365621
[17]:
from lapy.diffgeo import compute_geodesic_f
from lapy import heat

tria = T.boundary_tria()
bvert=np.unique(tria.t)

# get heat diffusion
u = heat.diffusion(T,bvert, m=1)

gu = compute_geodesic_f(T,u)

cutting = ['x<0.5']
plot_tet_mesh(T,vfunc=gu,plot_edges=None,plot_levels=True,cutting=cutting,edge_color='rgb(50,50,50)',html_output=False,flatshading=True)
Found 4800 triangles on boundary.
TetMesh with regular Laplace
Matrix Format now:  csc
Solver: spsolve (LU decomposition) ...
TetMesh with regular Laplace
Matrix Format now: csc
Solver: spsolve (LU decomposition) ...
Mesh is oriented, nothing to do
Found 3040 triangles on boundary.